#pragma once

#include "PLMatrixCoordinate.h"
#include "PLVec2.h"
/*#include <iostream>*/
#include <map>
#include <vector>

/**
  *	Stores the eigen vectors and eigen values of the passed matrix into the
  * passed std::vectors
  *
  * @param matrix matrix, for which the eigen vectors and eigen values are calculated
  * @param startVector start vector for the arnoldi iteration
  * @param n number of rows (and columns) of the matrix
  * @param eigenVectors pointer to an std::vectors, that is used to store the calculted eigen vectors into it
  * @param eigenValues pointer to an std::vectors, that is used to store the calculted eigen values into it
  * @return true
  */
template < class Matrix, class Vector >
bool GetEigenVectorsAndValues( const Matrix &matrix, const Vector &startVector, const PLuint &n, std::vector<PLVec2d> *eigenVectors, std::vector<double> *eigenValues )
{
	Vector *q = new Vector[n];
	
	q[0] = startVector;

	const Matrix &A = matrix;

	map<PLMatrixCoordinate, double, PLMatrixCoordinate::PLMatrixCoordinateCompare> H;

	int j, k;
	double value = 0.0;

// 	int count = 0;

	// Arnoldi Iteration
	for (k = 1; k < n; k++)
	{
		q[k] = A * q[k - 1];

		for (j = 0; j < k; j++)
		{
			value = dot(q[j], q[k]);

			if (value != 0.0)
			{
				H[PLMatrixCoordinate( j, k - 1 )] = value;
				q[k] = q[k] - H[PLMatrixCoordinate( j, k - 1 )] * q[j];
			}
		}

		value = norm(q[k]);

		if (value != 0)
		{
			H[PLMatrixCoordinate( k, k - 1 )] = norm(q[k]);
			q[k] = (1.0 / H[PLMatrixCoordinate( k, k - 1 )]) * q[k];
		}

// 		count++;
// 
// 		if (count == 10)
// 		{
// 			std::cout << "k: " << k << std::endl;
// 			count = 0;
// 		}
	}

	delete [] q;

	return true;
}

/**
  *	Returns true, if the matrix is positive definite, false otherwise
  *
  * @param _matrix matrix, for which the the definiteness is calculated
  * @return true, if the matrix is positive definite, false otherwise
  */
template < class Matrix >
bool IsMatrixPositiveDefinite( const Matrix &_matrix )
{
	int n = _matrix.dim( 0 );
	double val = 0.0;

	//Matrix matrix( _matrix );
	std::map<PLMatrixCoordinate, double, PLMatrixCoordinate::PLMatrixCoordinateCompare> matrix;

	// fill the matrix
	for (int row = 0; row < n; row++)
	{
		for (int col = 0; col < n; col++)
		{
			val = _matrix( row, col );

			if (val != 0.0)
				matrix[PLMatrixCoordinate( row, col )] = val;
		}
	}

	double sum = 0.0;

	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < i - 1; j++)
		{
			sum = matrix[PLMatrixCoordinate( i, j )];
			
			for (int k = 0; k < j-1; k++)
			{
				sum = sum - matrix[PLMatrixCoordinate( i, k )] * matrix[PLMatrixCoordinate( j, k )];
			}
				
			matrix[PLMatrixCoordinate( i, j )] = sum / matrix[PLMatrixCoordinate( j, j )];
		}

		sum = matrix[PLMatrixCoordinate( i, i )];

		for (int k = 0; k < i-1; k++)
		{
			sum = sum - matrix[PLMatrixCoordinate( i, k )] * matrix[PLMatrixCoordinate( i, k )];
		}
			
		if (sum <= 0.0)
			return false;				// A ist nicht positiv definit
		else
			matrix[PLMatrixCoordinate( i, i )] = sqrt( sum );	// Summe ist positiv
	}

	return true;
}